close all
clear
clc

% Common value for simulation
c0 = 299792458; % light speed
nmax = 1; % free space
freq = 1e9; % 1GHz
% LAMBDA = 0.08; % wavelength
LAMBDA = c0/freq;

Zlength = 5; % total Z for simulation

% =================================================================
% Start calculate grid resolution
NRES = 20; % resolve the wave with at least 10 cells, better >=10
NDRES = 4; % normally 1~4 for resolution of feature size
dmin = 0.2;

 % Compute default grid resolution
dz1 = min(LAMBDA)/nmax/NRES; 
dz2 = dmin/NDRES; % dmin means the smallest feature size
dz = min(dz1,dz2);
% Snap grid to critial dimensions
N = ceil(Zlength/dz);
dz = Zlength/N;

Nz = floor(Zlength/dz);
% End calculate grid resolution
% =================================================================

% =================================================================
% Start calculate delta_t and tau 
Nt = 10; % Nt >=10
tau = 0.5/freq; % tau <= 1/(pi*freq) 

% Delta T resolution
% dt1 = tau/Nt;
dt2 = dz/(2*c0);
% if using absorbing boundary, need make sure dt = dz/(2*c0),
% in order to make sure each cell need two time step

% dt = min(dt1,dt2);
dt = dt2;
% End calculate delta_t and tau 
% =================================================================

% =================================================================
% Start calculate simulation time
t0 = 6* tau; % at least 6 tau for a whole pulse
tprop = Nz*dz*nmax/c0;
T1 = 12 * tau; % total simulation time should >12 tau
T2 = 5 * tprop; % total simulation time should > 5 bounces
Ttotal = T1 + T2; % for highly resonant devices, this is not enough
% End calculate simulation time
STEPS = round(Ttotal/dt);
time_axis = (0:STEPS-1)*dt;
% =================================================================


% Initialize Materials to free space
ER = ones(1,Nz);
UR = ones(1,Nz);

nz1 = 100;
nz2 = 150;
ER(nz1:nz2) = 4;
UR(nz1:nz2) = 2;

% Compute updated coefficients
mEy = (c0*dt)./ER/dz;
mHx = (c0*dt)./UR/dz;

% Initialize Ey and Hx to zero
Ey = zeros(1,Nz); 
Hx = zeros(1,Nz);

E3=0; E2=0; E1=0;
H3=0; H2=0; H1=0;

% calculate source for SF/TF calculation
t=(0:STEPS-1).*dt;
Esrc = exp(-((t-t0)/tau).^2);
%Esrc=0.8*sin(2*pi*500e6*t); % 500Hz single tone for test
delta_t = dz/2/c0 + dt/2;
Hsrc = -exp(-((t-t0+delta_t)/tau).^2);
%Hsrc = -0.8*sin(2*pi*500e6*(t+delta_t)); % 500Hz single tone for test

% =================================================================
% Calculate Fourier Transfer
NFREQ = 1000; % seperate freq to 1000 point for FFT
freq_fft = linspace(0,freq,NFREQ);

% initialize fourier transforms
K = exp(-1i*2*pi*dt.*freq_fft);
REF_fft = zeros(1,NFREQ); % reflection
TRN_fft = zeros(1,NFREQ); % transmission
SRC_fft = zeros(1,NFREQ); % srouce

REF_t = zeros(1,STEPS);
TRN_t = zeros(1,STEPS);
SRC_t = zeros(1,STEPS);

% =================================================================

fig=figure;
set(fig,'Name', 'FDTD 1D SF/TF Simulation');
set(fig,'NumberTitle', 'off');

Nz_src = 50;

% Main FDTD Loop
for T = 1 : STEPS
    
    % Update H from E
   for nz = 1 : Nz-1
       Hx(nz) = Hx(nz) + mHx(nz)*(Ey(nz+1) - Ey(nz));
   end
   % Dirichlet Boundary Conditions
   % Hx(Nz) = Hx(Nz) + mHx(Nz)*(0 - Ey(Nz)); 
   % Absorting Boundary Conditions
   Hx(Nz) = Hx(Nz) + mHx(Nz)*(E3 - Ey(Nz)); 
   % update SF/TF in H field
   Hx(Nz_src-1)=Hx(Nz_src-1)-mHx(Nz_src-1)*Esrc(T);
   H3=H2; H2=H1; H1=Hx(1);

   % update E from H
   % Dirichlet Boundary Conditions
   % Ey(1) = Ey(1) + mEy(1)*(Hx(1) - 0); 
   % Absorting Boundary Conditions
   Ey(1) = Ey(1) + mEy(1)*(Hx(1) - H3); 
   for nz = 2 : Nz
       Ey(nz) = Ey(nz) + mEy(nz)*(Hx(nz) - Hx(nz-1));
   end
   % update SF/TF in E field
   Ey(Nz_src) = Ey(Nz_src)-mEy(Nz_src)*Hsrc(T);
   E3=E2; E2=E1; E1=Ey(Nz);

   
   % update fourier transforms
   for nf = 1:NFREQ
       REF_fft(nf) = REF_fft(nf) + (K(nf)^T)*Ey(1); % check reflection in k=1
       TRN_fft(nf) = TRN_fft(nf) + (K(nf)^T)*Ey(Nz); % check transmission in k=Nz
       SRC_fft(nf) = SRC_fft(nf) + (K(nf)^T)*Esrc(T); 
   end

   REF_t(T) = Ey(1);
   TRN_t(T) = Ey(Nz);
   SRC_t(T) = Esrc(T);
   
   
   %Plot
   imagesc(1:Nz,0,ER);
   % colormap('jet');
   hold on
   
   plot(Ey,'-b','LineWidth',2);
   hold on
   plot(Hx,'-r','LineWidth',2);

   title(sprintf('Step: %d of %d',T, STEPS));
   xlim([0 Nz])
   ylim([-1.5 1.5])
   drawnow;
   hold off
   %pause(0.01);
   
end

REF_fft = REF_fft * dt;
TRN_fft = TRN_fft * dt;
SRC_fft = SRC_fft * dt;


% -------------------------FFT Plot----------------------------
% Here is fft plot, however, some gaps and don't know why
%    - why need /dt? when calculte amplitude

sampleRate = 1/dt;
samplePerPeriod = sampleRate/freq;
signalTimeLength = STEPS * dt;
signalPointLength = signalTimeLength * sampleRate;

plot(freq_fft,abs(TRN_fft)/signalPointLength/dt*2)
plot(freq_fft,abs(REF_fft)/signalPointLength/dt*2)
% -------------------------End----------------------------



